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1 Introduction 

- -"irir;*": — - — - «» «~ 

-znrrr;*— • 

the P< “ «*>„ order most * increa J 
derivative approximations mus , be pe rfomied with Fas[ ' ' ^ ^ lhe 

efficient. If matrix mn,tip, icalion is used ^ ~ - * 

— of decrees of freedom to He practical “° ^ ~ * 

severe since the time step deer * 6 St6P restriction s are 

• - — - .. „ 

rr- — • - - — ... .. 

iuing the computational - 

subdomains, on which the spectral a • mt ° muUlpIe zones > called 

me spectral approximation is aDnlied Ac 

be used on more complex • P ‘ As a result - method can 

mpiex geometries. The use nf i„ 

polynomials in each subdnmfl- WCr ° rder a PP r °ximating 

7 “ 1 * “ 

:r;: -— «■ - - — 

Less than a decade after they were first im h 
multidomain methods for that hav h uced f24J, the bulk of the spectral 

h u h h b6en P r °Posed compressible fln 

hyperbolic problems still define the ,■ ° r similar 

for general hyperbolic problems, and [26J for the “ ' ‘ *' ^ ^ I29 ‘ ^ K 

Methods for the advecive terms of the compressible “ 

presented in [181 and TlQi a ■ v ier stokes equations were 

■~ns .as proZl;;™^ ^ ^ ~ “ — - 
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, hf> inhatto arid methods are whether the equal, ons 
The main differences between - ^ ^ „„ jn whlch the 

are written in conservattve or non-conserva ^ ^ was usedi (or eM mple, in 

interfaces are treated. The conservative form o cons idered in [331. 

, i ri«iandl41 Non-conservaUve forms 

,he methods presenmd the conserv ative form of the equations 

may lead to loss of conservation. interfaces to ensure that waves 

r— 

propagate properly through • norma l derivalive are 

*1 * m At least two values ui 

considered in more detai in - subdo mains that have that point 

available an interface point, & differenlial compatibility equation for the 

in common. One interface met o Derivatives are chosen from appropriate 

points along the interface , • inded ”. The other approach uses a 

subdomains so that wave To nl the correction method, the 

correction procedure «*1. MA ^ ^ at me boundaries. As a 

interior point approximation is integrated ry characteristic 

multiple solution values are available at e, nte^ P ^ ^ ^ , 

combina, ion of these so, u, ions is then made to correct 

waVeS aCr0SS 1,16 taBrfaCe ' ( has its adva „tages and disadvantages. Integrating the 
Each interface treatmen accuracy 

• that the solution can be approximated to any 

differential equation means th integration 

u - a rvf thp time integration scheme. impi^ 

in time, depending on the choice tangential space 

schemes can also be use U- tprfaC e S in more than one space 

derivatives must be 0 f geometnes on which so, noons 

dimension. This requirement sev ^ (he transf ormations of the 

can be computed, since „ means ^ be conlinuous across 

mappings between the subdomains and 
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subdomain interfaces p jn-_ 

smoothness of the grids ^ °" “* ° ll,er "“ d - do « "»> requtre 

However, the tempo,*, ^ ^ 

[7], page 245.) 0ITeCU “" me ‘ h0d 1S hmited '0 first order. (C.f. 

A disadvantage shared bv the two t„,. , 

method is simple to apply i„ 0 ' “* ' rea ‘ me,US * thdr com plexity. Either 

must be mad a, ^ ^ ~ ~ d — s a chotee 

— s most r wwch s — - — - 
approximations a, the comers of hd alg ° ri ‘ hmS Ca " * deVd ° ped f ° r ,he 

averaged grid origi„ a ,| y proposed by ^ °" ““ CMySh ^ 

quantities are defined on the Gauss-Cheh a ." melhod - “cell" averaged 

uaual Lobatto points [35], [15], (16 , ^ ^ Wl ”‘ e fl “ XeS *“ defmed at ,te m <™ 

disadvantages of the m e, ho ds Jus , descn J nTsTly me ' h0d ^ ma " y * ' he 

approximated to any temporal order of accuracy „ h , a "« it can be 

il do “ not requtre conttnuity of the transf * e ““ e Wcall y flexible because 

y ne transformations across inrprfa^c t 

apace d, mens, on, the method does requite special attention a, the ‘ " m ° re ‘ ha " °" 

Currently, a simple average of the at the. comers of subdomatns. 

contributing subdomains [16], ^ ^ S °' UUOnS “ C ° mpU,ed and broadcast to all 

■u this paper, we present a new multidomain spectral col, • 
solution of compressible flow problems The co "oca„o„ method for the 

analogous to fully sta d n,elh0d “ based on a staggered grid, 

y ggered 8nds often used with finite 

soludons are defined at the nodes of a r a meth ° dS ' The 

^ Gauss quadrature rule and n 

nodes of a Gauss-Lobatto rule. Staggered-.rid ’ n “ XeS are evala ated 

proposed for the solution of the ' “ spectre! approximations were first 

solution of the incompressible Navier-Stote, , 

234.) Our grid will be tdentical to the f „ " S ' <C f ' (7) ' P a * e 

-ucalto dte fully staggered g ri d of Bernard, and Maday [2], 
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for compressible now problems has all 
The staggered grid mul.idomam metho ^ ^ avera ged 

f «ri in the methods discussed above, r\ , 
the desirable features 10 ^ be possl b,e to a PP , shock capturing 

method it is conservative. • d ind ependently ot their 

techniques to the approximation. Subdoma ^ inKrfaC e condition can be 

tHp method is geometrically fl 

neighbors, so the method However, in multiple space 

imputed to the same temporal « - * " ^ ^ _ of 

dimensions the method does not me u corn ers, and . 

any number of subdomams can meet at ^ in next section for 

The paper, divided a ^ of notation used throughout 

problems in one space dimens , conser vative, while methods that 

t W r. Wesimw and a, i^rsysmm Will, used as 

upwind derivatives are . . ^ ^ for ^ problems, 

examples to show problems, an example is 

- - - — - 1 — 1 ,■ »- v 

u Vs-rvh order temporal accuracy can also be oDtame 
included to show that tag ^ ^ ^at the method remains 

« describe Ihe algorrrhm in rwo space dim ^ ^ ^ of the 

— - » - l . - - ■ - 

use of the method for two-drmensro ^ ^ ^ eiponentia i accuracy is 

source flow, for which there is an exact so a circular bump 

K , m The second problem is a subsonic now u 
Obtained for this pro ^ ^ ^ ^ decay exponentially fast Finally, we 

i„ a channel, and we s o converging . diverg , ng norrle and compare die 

solve a transonic flo de in ^e Section 5. 

results to experimental data. Concluding remarks ar 
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2 JEST" G ' id APP,0Xima,i °" "> SPaee 

2.1 Notation 

The staggered grid approximation uses two gride 
and advective Huxes Uni t h compute the solution values 

only the nodes of th G C ° mm0n ™ uses 

- Mh die Gauss and ^ ~ — 

^USS-LobattO points. We dpnnto 

by the Lobatto points X and the r - P 0l nts on the two grids 

puinis, Aj, and the Gauss points X 

’ >+1/2 > 


X ‘ 2 f 1 C<JS f 77 j j j = 0.1 N 

w£[i-w(iZ±L,)] 


<D 


I" OX we have mapped the usual collocation points defined 

convenient unit interval. The overbar and half I • J to the more 

used only for their value P0 ‘"‘ n0ta ‘ i0nS ** ‘ he GaUSS ^ ™ 

niust be understood I “ ^ ^ ~ ,t 

[71 . 

Two polynomial approximations are defined one f „ 

polynomials of degree less than nr i F eacl1 ® nd Lel ‘he space of 

* than or e q ual to IV be denoted P„ Le, be , he 

Lagrange interpolating polynomial 

',(i)= n| 

111 v - ' (2a) 


i=nfizi 


1=0 l Xj — X ; 
**j 


the Lobatto grid. On the Gauss grid, we define h eP t0 „ „ 

polynomial ‘" n e P „_, t0 be the 


N-if 


7 + 1/2 


i«)=n 


S-x, 


f=0 l * 


7 + 1/2 


y+1/2 X i+m 


(2b) 
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M , - ““ - - 5 “ “ ‘ * 

„ ^ ™ - - « “i— " — 1 ■ 

n (3a) 

<2(X) = 2M (X) 


=0 

jV-1 


Q(X) = Y J Q i+ u2 h j + inS X> >- 

;=° 


(3b) 


llThvOM.Dimrmlima! Sum"* ‘‘" ,l ' w Stator Ef«c«i> 

T » ** ~ - •“*' " 

scalar problems of the form 

(u, + />) = » 9/ / > o, x e [a.H (4) 

■ m(x, 0 ) = Mq(-x) 
u(a,t ) = g(0 

« M, ^ 

t _ ! 2 * wh,ch are ordered left to right. A simple linear dansfomtauon can 

l0 a* „„U interval, so that on each subdomain we solve the probiem 

u , + i-/» = 0 Xe [0.11.O0 (5> 

x x 

el „rid defined bv (D- For convenience. 
On each subdomatn is placed the staggered gnd defined y< 

------ ™;:r :r:ri:— - 

reauired by the method. We then let V 

exact solution, u on fl». Similarly, the flux is approximated by the P<r yn 
F* (X) e P w , defined by (3a). Substitution of these approx, mat. ons mto 


Tjk \_dF\X)_ = /? t( X ) k = \,2,...,K 

U l ^ „ 3Y 


( 6 ) 


x. ax 
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To obtain the equations that define the solution unknowns at the G 

that the residual, R, be zer o at the r ° P ° imS ’ We require 

collocation app^J " ^ ^ “ — * -ds to the 


dU 


7 + 1/2 


+ 
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dt x x dX J "’"’-^v-i ( 7 ) 

The spatial derivative operation in n\ 

'*** of flux values by a derive matrix. D 1^“ “ ^ mUmPiiCa ‘ i0 " ^ 

dF k (X ) 

'7+1/2/ 


_ 7+1/2/ . n 

Km )K = Yd I ,F 

n=zQ / J n r 

n =0 


( 8 ) 


so we write 


dF k 

dX 


0 + 1/2 


: ( DF ‘U=I«. 

*=o 


( 9 ) 


Thus, (7) can be written in vector form as 

d U* 


dt 


+ DF k =0 k = 1,2,... 


GO) 


where U* = \rj k u k n* i r + r 

l 1/2 U w - U n- U2 ] ,F* =[F k F k ... F k J r , 

To compute the flux values on the r k 

reconstruction procedure We f ° aU ° gHd ’ WC USe the lowing 

p ocedure. We first evaluate the interpolant U k (X) e P at thp . h 

P° lntS by muW ^ «* vector of solution values „ y m ’ "" L ° baU ° 

es by an inteipolation matrix, I. 


N~\ 


U (X ; ) = U u ( y \ V 

J 4 L/ n+l/2 n n +U2\X t ) = > /. TJ 

*=0 ; 4-1 J,n+l(2 U rt+ll2 (H) 

The family of characerist.es of (5) runs left , 0 ri » ht " l„ s 

— - - 'eft subdomain boundary!, ^ J * “ * * 

«— t “ 
to define the j = 0 value on the r u h boundary condition 

two values ^ (1) ^ 7 “ S “ bd ° ma,n ' * ^ where 
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■ f e The result is an upwind evaluated approximation at both the left 
boundary and the interfaces. The fluxes, b,, 

on the Lobatto grid. through the definition of 

The method imposes the boundary conditions, weakly, through 

The metnoo v at the boundary or 

the flux, since the discrete solution values are ^ ^ 

. i t u p cin^le domain approximation of ( ) J 

interfaces. ' lerras „f the interpolant UW and the 

we can write the flux F(.X)-U » 

boundary condition as 

_ ( 12 ) 

U{X) = U(X) + [g- v ( x > ’ 


so that the polynomial U(X) satisfies 


UiXj) ■ 


U(Xj) ; = 1,2,...,N 

U i=° 


(13) 


Then (7) can be written as 

+ -^-U'{X j+u2 ) = ~~ [U\X 0 )- gVo(Xj+m ) 7 

’—■“"“T. , „„ „ - “ 

Equation (10) is a system T nnrinciole 

, „ in time to get the approximate solution values at die Gauss points. In prmcip , 
integrated in time to g We have chose „ , 0 use low storage 

any common integration proce ure c ^ ^ for ^ computation of 

Runge-Kutta methods that require is on l v an iterative procedure, we 

steady-state problems, for which the time discreuzauon ,s only an 
have used a mid-point rule, 


77*,rr+l/2 _ T/<=." 

U l+U 2 ~ u i+ 1« 


A f _ 1 _ dF k ' n 

' 2 x Y ax 


lj+1/2 


rri.n+i 

U j+U2 


■ TJ k ' n 

- U j + 1/2 


-At 


1 dF 


.lt,n+l/2| 


ax 


j = 0,1 N-l 

; = 0,1 N-l 


(15) 


b+i/2 
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temporal"!! '° haVe * ^ ^ Slep lhat can used ^ 

emporal damping rntroduced by ^ scheme ^ ^ ^ 

eigenvalue structure of the differentiation matrrces, other chorees m, g ht Include schemes 
optimized lor raprd converge to s, estate, such as those drscussed in (I0) and 

For time dependent problems, we have used the third order 2-N storage method of 
■ tamson 9], and the more recent fourth order scheme of Carpenter and Kennedy [8] 
note that .. should also be advantageous use the new low storage Runee-Kutta 

methods derived by Hu „ * 121,. wh.ch are optimized to min.mize the phrase and 
drss.patron error rntroduced by the temporal approximation. 

for the ^alar problem deam^ed abovef ^ Pr ° CedUre ' ^ ^ a, * onthm 

Algorithm I. (staggered Grid, Scalar , 1D) 

1. Interpolate u = fu u 77 l r 

L 1 / 2 ’ 3/2.-.U n _ 1/2 J to the Lobatto points: 

Compute the matrix-vector product U* = I*u* defined hv m w u 

v aerined by (1 1) for each subdomain 

^.Compute the flux vainee *. ■ 

F k ~f(U k ) ■ 3 lnternal P°ints on Lobatto Grid- 

| 3 . App ly ^boundary and interface conditions: 

F o=f(Kr ] ) k = 2,3,... ,K 

4 - Di — “* — - « O— -id by subdomain: 

mpute the matrix-vector product D*F* k = 1 .2 N by eq. (9). 

|S. Update the solution by subdomain: 

Integrate (.0) by the chosen ODE solver, repeating Steps 1-4, as necessary. 

|6. Repeat 1-5 until done. 
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we note that the method rentes two matnx-vector products pet Runge-Kutta 
stag , This is twice the wo* of a Lohatto grid method, ot of the ceii-avetaged method. 
Thus, mere is no speed advantage for the method in one space d.menSHnn ^ ^ ^ 
A desirable feature of the staggered gnd approximation, ), 
conservative. To show conservation, we define the quadrature 

1 N — 1 

f F(X)dX = £ F j+ mWj+m VF G P *-' 

o >=° 

w ; - +1 ,=i w**® 

. ry. sum over all points and all subdomains is 

For each j, we multiply (7) by x x w j+v2 - 


(16) 


Jpx-J f t dUj+m , p'k 1 = 0. 

dt i+m ) 

1=1 i =0 ' 


(17) 


Now, 0*(«.f-(X)X. and x x is a polynomial of degree zero, so we can replace the 
sum over j by integrals to get 


il 

1=1 o' 


^^+F*(X)V = 0 
dt J 


(18) 


Upon integrating the flux derivattves, the tuterface contributions cancel, and 

iI|j[/‘(X,x*<«| = F ,( 0)-F' 1( l) (,9> 

In contrast, a mlidomatn method defined on the Lohatto points, where the 

upwind value of me flux — * - - - “ ^ “® h ^ 
■ To show this it is sufficient to consider two subdomains, tf- and Q 

system of equations 
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j = 1,2,... ,N -l 


U o =g 
dU_ L 

dt 

dU R 

dt =~ F ; R j = l,2,...,N-l 

dt dt ~~ Fn 
dU* 


~ -f' L 

j 

1— - —f' r 


( 20 ) 


dt 


iL — __ rp f R 
r N 


We then multiply each equation by its associated Clenshaw Curtis • u 

get ^enshaw-Curtis weight [7] and sum to 


f du } L L ^dUf . 
j=o at ^ dt j 


+ </J % w ‘ F ‘ L '% w i F ’i R ~<[F-„ L -K'‘] 

= ^(0)-F«(l)- % . [r/ _ F ., ]+ y L + ^ 

V dt 


( 21 ) 


~ - - - - - 

problem stnce the dlfference be(ween (he “ " 0t expec “ •« * a 

fas,, while Ihe coemcien, decays as 0(l/*2* 

As an example of solutions computed asms Algorithm I u, 
solution of the equation u +u - f r e ,0 ?, „ *' ™ hm 1 we c ™P«e a steady 

a ' ‘ °.21t>0.Scalarexamplesofo„edimensional 

time dependent problems can be found in [271 m 

difference methods can be found [37] ‘ “ “ * "** ° f — 

chosen so that the exact steady solnUon is ^(<7- "" 
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Fig. 1. Steady solution of a scalar wave equation 

Convergence of the error rs exponent for this problem- B, 2 shows fhee.or 

• i A&r fnr the subdivision shown in Fig. * 
pl0 „ed as a function of the Gauss polynornta or raelho d 

For contparison, we have aiso plotted the error of the srng ^ _ 

(20 ). We see that the staggered grid appro— , a. ieast 

staggered grid approximation, and sometimes more accurate by a factor o 
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yStemS 0f hy P erboJic equations of the form 

I i IT" / y\ v „ 


[Q< + K(Q )-0 X(=[ a ,b],t>Q 
(QU,0) = Q oW 


( 22 ) 


where Q and F are w-vectorc u, 

Jacob.au matrix A =a F/aQ ^ ^ ^ * hyperbolic ’ that is, the 

assume, as is the case tor the L ’ ^ & ^ ^ We ^rther 

for the Euier gas-dynamics equations, that the flux can be written 
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conditions are applied. equation, except for 

nf the system follows that of the scalar equ 

The approximation of M an interface between subdomains 

rrr rrrr;— — «■ - * 

r::- — — — — — 

F - AO = ZAZ'Q = ZA'Z 'Q + ZA Z ‘ Q (2 ’ 

- a + ,M The first term represents waves moving left to right, and the second 
where ^ Anupw i„dapp—n chooses flf - the 

represents waves movtng . compo „enls to give 

right going components, and Q„ for the te % 

•. _ . rw k~ r 7~ l r\ k (24) 


Fr = n = 7(q‘» , .q‘.)* za * z " q » ,+za " z " ,q ‘°' 

extensively in the finite dtfference commnntty (Ref. PO). 

vector splitting and flux difference splittmg. ^ed using 

The resolutron of the jump at the subdomarn mmrfaces 

, „ The nux is decomposed into a right gorng and a left g 
flux vector splitting. The fl 0 f F+ has only 

_ + F- (O). The splitting is done so that the Jacom 

F(Q) (Q h - of F - has only negative eigenvalues. Examples are 

positive eigenvalues and the Jaco ian - Steffen [30]. Using 

negative flux is computed usrng the value from t 

^.F!-?(QSr‘.<U)-^) +r ( Q 5) ' 

a ■„ ti«l to compute the flux derivatrves at 

Van Leefs Flux vector spUtting was used » U« 

art nf thP Navier-Stokes equations, 
interfaces for the advecuve part of the Navte 
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-:“;rrr *• - — — - - 

—» m ~r — — - *. 

fluxes. However, we found tha , the Va „ ’ °" ly * S ‘ mP ' e S ™ °' ‘ he 

^enre was unstable for some long Ume integltiol ^ **“ 

~ ;::r * . . 

KamiadakisflS], Several solver choices! ^ a ™"* 0aI0S a " d 

with the entropy fix Formally o‘ u C ave use( * fl° c ’s [34| solver 

™ aJly ’ 81Ve " ,he «“» Q»‘ and QJ. we wrile 

^ Q * ’ Q ^ = l( F ( Q ^) + F (Q:))-|R|A|R-'(Qj_ Q r-.j p6) 

where R is the matrix of the right eigenvectors of the Jacobian of F 

Rop-average of Q k ~ l and o* Th F ’ com P uted using the 

points^ Qo- The eq. ( 26 ) is modified to correct the entropy across sonic 

the solution assumed k! * TT " ^ C ° mPUted “"*» a "« 

- — ^-—--.-.huown,^. 

F l = r(Q(a.», QJ) 

on the left, and 

F i - rfQl.Q (t,D) 

on the right, where Qfa, t) and Q (h t) re . re (2?b) 

Other ways to compute the h a eXtenOT S0 ’ Uti ° n * ‘ he bou,Klll ™s. 

y u com Pute the boundary flux wh^n n, a f 

will be described in regard to specific probl ■ ” 6X16001 S ° 1Uti ° n “ ^ kn ° Wn 

s ° specinc Problems in Section 4 

1,1 SUmmary ' f ° r SyStemS ° f « bave - following algorithm: 
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Algorithm II. (Staggered Grid, system, ID) 

Isolate the Gauss-point sedation values to the hohatto 

Compute the matrix-vector product Q* = l‘Q* ^ eq. (H) for each subdomam. 

| 2 . compute the interior point £lu*es: 

F*=F(Q‘) ) = l,2,....N,t = 1.2 K 

| 3 . App ly the interface conditions: 

f‘„- , =f‘ 0 = 5(Q» , .<IJ) k=2 -- K 

L . App ly boundary conditions at left and right 

■ i derivatives at Gauss points: 

1 4. Compute spatial d 

Compute the matrix-vector product by eq.(). . 

5. update the solution at the Gauss points 


d Q 


kl + F'‘,„=0 j=0X...N-l,k = l,2 K. 

dt 1 


repeating Steps 1-4 for each Runge-Kutta stage. 

[6. Repeat Steps 1-5 until don 

As an exampie of die application of Algorithm II, we solve the problem 
Q i + F,=0 x e[-l,4], t > 0 


(28) 


Q = 


M F-f 1 2 lc 

v’ F ' 2 IP 

L J 


Further examples that mode, unsteady acoustic propagat.o 
propagation in a quasi-one-dimensional nozzle, can be fou ™* ^ 

conditton for (28) was chosen to be a Gaussian pulse with the peah at x = 

f -12 (x-1) 3 1 


(29) 


Q(*>°)= o 


we specify boundary coupons so that the waves pass through the boundanes without 
reflection. 
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- v(-lt) = e -> 2 ('- 2 ) 2 
u(4,t) + v(4,t) = e -^(3-3„ 2 


( 30 ) 


same number ° f w ithin l su zr nr ^ and with ihe 

- - - * — , as — — - y 

computed the solution on the finest grid A t - „ . StU<ly ' he ,emPOra ' we 

rounding error. A plot of the etrors as a f ‘ “ ‘ he ^ was c 'ose •» 

the th,rd ° rdW «» for the fourth order tnedrod ole ^ ‘ ^ " 2 " 5 “* 
accuracy ,s obtained. We aiso see tha, over the A 1 “ ^ ^ 

observations of [8], ' ThlS 1S COnsist ent with the 
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log (Maximum Error) 



Fig. 4. Decay of the spatial error of the system (28). 





3 The Two-dimensional Approximation. 

We now describe the approximation 
conservative form, 


Of the Euler equations of gas-dynamics 


in 


dQ dF dG 
dt dx + By ~°’ 


(31a) 


where Q is the vector of solution unknowns 


vectors 


and F(Q) and G(Q) are the advective flux 


Q = 


P 

pu 

pv 

pe 


pu 


G = 


pv 

puv 

p + pv 1 
v(pe + p) 


(31b) 


F= P+pu~ 

puv 

u(pe + p)\ 

We assume Y = 1 .4 and that p e = p / ( y _ 1} + ,2 2 , 

such th ( V )/2. For axisymmetric problems, 

e ransomc flow in the converging-diverging nozzle discussed later we 

T! ' he ““ C00rdi " ate a " d '' " ' he rad,a ' — We .hen add th e tight 
hand side of (31a) the vector n * lu 


H = 1 

}’ 


pv 

puv 

pv 2 

v(pe + p ) 


(32) 


3.1 Mapping in two space dimensions. 

In .wo space dimensions, we subdiv.de , computationa] ^ Q _ 
quadrilateral subdomains. il l - k — \ 2 k f- 

region inm f „ ^ 6 Shows ™ * a division ofa 

egion into four subdomains We mate .a.,.., 

paper Firs, „ P,,0 " S ab0UI “* Vision <« 

paper, first, we allow subdomains lo intersecl onlv , ■ 

Secon , , °" ly at a P° ,nt or along an entire side. 

Second, we assume that the approximation is conforming so that grid r 

across subdomain interfaces Finn,, C °' nC ‘ de 

"terraces, finally, we assume that the subdomain boundar.es do not 

P 7 ■" tirae ' - — we Will mahe the assumpiion that the II 

“ ° rd “ ‘ S “ Sed " - - - * each subdomain, in practice! 
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number grid points can vary, as 


long as the approximation remains conforming at 



Fig 6. Diagram of a subdomain decomposition in 2D 


Subdomains are mapped onto Ore unit square by an isoparametric mapping. Let 
the vector function g(s), 0< t Si define a parametric curve. Define also the polynom 
degree N drat interpolates , at the Gauss-Lobatto points to be 

ns)=2,s( s i) < i (s) ' 13 ' 

counted counter-clockwise, bound each 
the unit square by the linear blending formula 

(34) 


;=0 


Four such polynomial curves, r m (s), m - 1,2,3,4, 

subdomain. We map each subdomain onto 

x N (x, Y) = (i - Y)T x {X) + it 3 (X) + a - x>r 4 (n + xr, (X) . 


_ Xi (l _ X)(l -Y)- x 2 X(l -Y)- X 3 xy 

of the subdomain, counted counter- 


where 


the x/’s represent the locations of the comers 


cW*isc. (34 ), the Euler equations (31) 

Under the mapping Cl K <-> |0,ljxiu,rj g 


become 


dQ + l 
dt J 


'djL + 

dG 

ax 

dY 

j 


= 0 


(35a) 


where 
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(35b) 


F = ^F-^G G = -y*F + x%G 

Sfnce we assume here that the subdomain boundar.es do no, move in time, we can write 
(35a) as 


^Q + ^F(Q) + aG(Q) 
dt dx dY 


(36) 


where Q - JQ and the fluxes are still defined by (35b). 

3.2. The staggered grid 

A fully staggered grid is used in two space dimensions. A schematic of the grid 

on a s *ngle subdomain is shown in Fig. 7. The grid is the same as the staggered grid 

proposed by Bernard, and Maday [2] for the solutton of the incompressible Navtet-Stokes 

equations. In what tollows, we will ignore superscripts that denote which subdomain is 
being considered, unless necessary. 

Y 




(i+l,j+l/2) 


Fig. 7. Diagram of the fully staggered grid in two space dimensions. 

Potnts of type “a” in Fig. 7 represent the Gauss/Gauss points 
o,o. M n)-‘-J-OX...,N-\. The g rid thal resul[s from these ^ ^ ^ 

product of the one dimensional Gauss grid defined in m um 

gnu aenned in (1). We approximate the solution 
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and the transformation Jacobian a, the Gauss/Gauss points, and denote them by 
_ . 7 - j(Y Y From these, we compute the Gauss 

Q/ + i/2i + i/2 i+1/: 2 ’ 1+U2 

_ 7 75 Finally, the interpolant of the solution 

point values Q i+1/2 ,; + i/2 - Ji+u 2 ,j+m'*i+u 2 . j +u 2 - 

through the Gauss points is a polynomial in P w -i.v-i = p n-i 0 P *~' \ 


Q (X,f)=XiQw«,j + W2^ 1/ 2 (X ^ +1/2(r) 

•a 


(37) 


i=0 j—0 


, t M : n Fi p 7 form the Lobatto/Gauss grid 

The points of type b in tig. / 

,y y w / - o 1 /v. On this grid are evaluated the horizontal flux vector, F an 

the metric terms yy and ay . The metric terms are computed as 3y ( „ j.m) 
d VfX,.^,,,,)/ 3T- At points interior to a subdomain, the horizontal flux is computed 

by 

where Q is a polynomial of the type (37) that passes through the values 

The computatIon of the fl “ at " y a " d interlace P0 ' nlS ,S 
described in the next sub-section. 

The vertical flux and the derivatives yx and ax are computed on the 
,, hv “r” on Fie 7. The points on this grid are 

Gauss/Lobatto grid, marked by on g. 

- V s . • A 1 N - \ The metric terms are computed as dy (X 1+1/2 , ,) 

dx s ( ^ +]/2 jY] )/dX. The vertical flux is computed at interior points by 




G w «j - >’* 


and at boundary points as described in Section 3.3. 

„ may appear that to define quantities on three different grids would lead to a 

complicated method than a single grid Lobatto approximation, Thts 

of the fluxes on the staggered grid by (38) 


significantly more 


turns 


out not to be the case. The definitions 


and (39) mean that the 


reconstruction procedure, i.e., the interpolation of the solution 


22 


needed to compute the fluxes at the Lobatto points, is not a two-dimensional operation. 
Rather, it is the less expensive sequence of one-dimensional interpolations, given by (11). 
The values of the solution vector required to compute the flux vectors are actually 


N-\ jV-L 


q<*,w=2;Iq 

i + 1 / 2 ,; + 1 / 2 ^/+ 1/2 )^;+l /2 (^y+ 1/2 ) 


i = 0 ;=0 
N - 1 


=Zq 


(40a) 


i= 0 


r+1/2, y+1/2 


2 W*,) 


and 




Q(-^i+l/2 » Yj ) ^ Qi+l/2,;+l/2^i+l/2 (-^r+1/2 )^; + l/2 ( ) 


i=0 j=0 

N-l 


S Qi+l/2.7+l/2^;+l/2 (^7 ) 


>=0 


(40b) 


since, by construction, 


^m+l/2 (^7+1/2) ~ ^m.j 

K+l/2 (-^i+1/2 ) = 


3.3 Interface and boundary treatment. 

To describe how we compute the interface and boundary conditions using the 
staggered grid approximation, we will refer to Fig. 8, which schematically represents four 
subdomains and the locations at which solution and flux values are computed. Only the 
collocation points near the boundaries are marked. The circles represent the solution 
values, which are located on the Gauss/Gauss grid. The locations of the horizontal flux 
values, F,. y+]/2 , are represented by solid squares. The locations of the vertical flux 
values, G i+1/2 ^, are marked by hollow squares. From the diagram, we see that along the 
interfaces between subdomains 1 and 2 and between subdomains 3 and 4, only the 
horizontal fluxes need to be computed. Along horizontal interfaces, like those between 
subdomains 1 and 3, only the vertical flux needs to be computed. Because the grid is 
fully staggered, the coupling is through subdomain faces only, not through the corners. 
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Fig. 8. Diagram of four subdomains showing locations near interfaces where solutions 
and fluxes are computed. Symbols: • solution; ■ ,F, □ ,G 

Fig. 8 indicates a significant advantage ot the fully staggered grid over an 
unstaggered grid. In the unstaggered approximation, for example as described in [26], 
special comer algorithms must be devised to ensure correct propagation of waves through 
the corners. Each special case must be coded separately. Also, the choice of bi- 
characteristics that determines the domains of dependence becomes more complex as the 
number of subdomains/boundaries that come together at a point increases, making the 
derivation of these special cases more difficult. The staggered approximation does not 
include subdomain comers, so conditions do not have to be specified at comer points. 
Any number of subdomains can come together at a point without the need for special 
point approximations. No special code is required even for very complex subdomain 
topologies. 

The interpolation of the solution by (40) produces two solution values at an 
interface point, one from each of the two contributing subdomains. As in the one 
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dimensional case, we do not expect these two values to coincide, except in the limit of 
infinite resolution. A single flux is calculated, as described for the one-dimensional 
problem, except that we only consider waves propagating normal to the interface. This 
normal wave approximation is common for finite difference approximations [20] and has 
been used in [18], [16] and [4] for spectral approximations. We note, however, that other 
two-dimensional wave decompositions are possible, like those surveyed in [32], 

Physical boundaries can be viewed as interfaces between the external flow and the 
computational region. Wall boundaries can be computed by imposing an opposing flow 
that enforces zero normal momentum flux across the interface. Subsonic inflow and 
outflow boundaries can be computed by replacing the solution that would have come 
from a neighboring subdomain by the free-stream values, if they are known. If the full 
state of the exterior flow is not known, the known quantities can be specified and the 
remaining quantities can be computed by a characteristic method. Once all solution 
quantities are known on the boundary, the flux can be computed. An example of this 
approach is provided in Section 4.3. Supersonic outflow boundaries require no extra 
conditions. 

3.4 Discretization of the equations. 

Once the fluxes are computed, the spatial discretization can be made. From the 
discrete flux values are defined the polynomials 

i —0 ;=0 

N-\ A f (41) 

G(X,Y) = XSG 1+ia ^, /2 (X)/ ; (f) 

i -0 j = 0 

Derivatives of the interpolating polynomials are then evaluated at the Gauss/Gauss grid 
points. Like the reconstruction procedure, the differentiation of (41) can also be done as 
a sequence ol one-dimensional operations 


25 



li+l/2,;+l/2 


(42) 


3F 

ax 


= £F.j + uiKiX^n) 


n- 0 


dG 
a Y 


i+\/2J+\/2 


=ld 




+1/2, V* ; + l/2 


n=0 


Because both interpolation and differentiation operations must be performed at 
each step, the total work of the staggered grid method is twice that of a method that only 
uses the Lobatto grid. The new method requires the same amount of work, however, as 
the cell averaged method [16]. The reconstruction procedure in two space dimensions for 
the cell averaged method is more complex than in one, and requires the same amount of 
work as both the interpolation and differentiation operations here. 

Finally, from the definitions (37)-(42), the semi-discrete approximation for the 

solution unknowns can be written as 


dQ 

dt 


i+m,j+U2 


+ 


3F dG 
dX + dY 


= 0, 

i+U2,j+l/2 


Eq. 43 can be integrated in time as described in Section 


i = 0,l,...,N-l 
7=0,1,..., N-V 

2 . 2 . 


(43) 


3.5 Properties of the staggered grid approximation. 

The staggered grid approximation is both conservative and free-stream preserving. 
A net gain or loss of Q is determined only by the flux through the exterior boundaries. If 
the solution is constant in space, then the solution must remain constant in time, even in 

the presence of a spatially varying mapping. 

We first show that the staggered grid approximation is conservative. It is 
sufficient to consider the four subdomains shown in Fig. 8. Let the quadrature weights 
w i + ii 2 ’ r lj+v 2 be defined so that 

^ , i+in,j+u2 w i+v2 r lj+in V/>€ (44) 

o 0 J~° 
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By exactness of the quadrature, the sum of eq. (43) times w, +l/2 7 ] y+1/2 over all the points 
within a subdomain is 


i i 


0 0 


^-dxdY= 
dt dt 


w 


(+ 1/2 


fy+1/2 


1+1/2, y + 1/2 


A/-1 


3F 3G 
+ 


= y 

ay 

1 1 


w, 


(+ 1/2 


fy+l/2 


(45) 


(+l/2,;'+l/2 


0 0 


9X 3F 


Thus, for each subdomain. 


d 1 1 i i 

— j J QdXJT = - J F(l, Y)dY + J F(0, Y)dY 

^ 0 0 0 0 

1 1 
-Jg(x,i)jx+ Jg(x,o)</x 


(46) 


When (46) is summed over all subdomains, the interior integrals cancel so that only the 
boundary contributions remain: 


d 4 l - 1 - 


— X JJ Q kdxdY = j (f 1 (0, Y) + F 3 (0, Y))dY - j (f 2 (1, Y) + F 4 (1, Y)}dY 


(47) 


J(G 3 (X,0) + G 4 (X,0))fi?X- j(G 1 (X,l) + G 2 (X,l))jX 


The staggered grid approximation is also free-stream preserving, which means 
that the isoparametric spatial mappings do not introduce false source terms. It is 
sufficient to consider the approximation within one subdomain, since all derivatives are 
computed locally by subdomain. If we take F(Q) = G(Q) = 1, then the approximation 
(43) becomes 


dQ 

dt 


+ 1 


(+l/2,;+l/2 


^ ( N N N \ 

(,>'}' + yx+ x x) 


dX 


= 0 


(48) 


J/+1/2J+1/2 


Since x e P 


N,N’ 
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(49) 


so that. 


d_ 

dX 


dx 

~dV 


N\ 


I+1/2J+1/2 


= f J xMx MI7 )r l (Y ltm ) 


3 

dr 


3x"_) 

y 1+1/2, j+1/2 ■ 


dx 


d Q 

dt 


= 0 


i+1/2 J+l/2 


1=0,1 A^-l 

; =0,1,..., W -1 


(50) 


4. Examples 

In this section, we use the staggered grid approximation to compute three steady 
flow problems. The first problem is subsonic flow from a point source, which has an 
exact, analytic solution. We use this solution to show that exponential convergence is 
obtained. The second problem is a subsonic flow over a circular bump in a channel. 
Although there is no exact solution for this problem, we show that the errors due to 
entropy generated along the curved wall decay exponentially fast. The final problem 
computes a transonic flow in an axisymmetric converging-diverging nozzle. That 
solution is compared to experimental data. 


4.1 Subsonic Point Source Flow 

As our first example, we consider the flow of a steady, irrotational flow exiting 
from a point. This flow can be solved exactly by a hodograph transformation [12], The 
streamlines are radial, and level curves of the Mach number, pressure and density are 
circles centered on the source. We will compute this flow in two geometries. The first 
represents a flow in an expanding duct, where two streamlines are chosen as walls of the 
duct. The second geometry, a square with five circles cut out of its interior, is included to 
show that the method can be used to compute a flow in a complex, multiply connected 
geometry. 

The first geometry represents steady flow in an expanding two-dimensional duct 
with straight walls. The lower wall was chosen to be the line y = 0 and the upper wall was 
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the line y = x tan(rc/6), for x between 1 and 1.5. The exact solution chosen sets the Mach 
number at the lower left corner of the domain to be M = 0.6. 

We examine solutions for three subdomain decompositions, each having four 
subdomains. Fig. 9 shows the three decompositions. In the first (Grid I), the subdomain 
boundaries are straight lines so that the mappings defined by (34) become bi-linear 
transformations. The second and third decompositions are included to study the effect of 
curved subdomains. Both perturb the Grid I by a sine wave of amplitude 0.1 into 
“bulging” (Grid II) and “wedging” (Grid III) decompositions, named so in [36]. 



Fig. 9. Three subdomain decompositions for the diverging duct problem. 

Wall conditions, applied as described in the previous section, are specified on the 
top and bottom boundaries. The left boundary is a subsonic inflow boundary. For that, 
we specify the exact solution as the incoming condition for the Riemann solver. The 
right boundary is a subsonic outflow boundary, and again the exact solution is used to 
specify the external flow. A perturbation of the exact solution was used as the initial 
condition. 

Fig. 10 shows the computed, steady Mach number contours for Grid I. In that 
figure, and in those following, contour lines are plotted using solution values interpolated 
from the Gauss points to the Lobatto points using (37). The grids in Fig. 9 show those 
Lobatto points. The solutions are represented interior to each “cell” bounded by the grid 
lines. The interpolation is done for display reasons, since a plot using the Gauss points 
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would show gaps between the subdomains, a result of the fact that the solution is not 
defined at the interfaces. Plotting the interpolant does give some visual indication of the 
size of the solution jumps at the interfaces. 



Fig. 10. Mach contours for flow in a diverging duct. 


The staggered grid approximation is exponentially convergent for the point source 
flow in the duct. For Grids I-III, Fig. 1 1 shows the weighted L 2 error in the density as a 
function of the polynomial order used in each subdomain. The most marked observation 
is that for this problem, the error and the convergence rate are not sensitive to the 
presence of curved interfaces. In fact, the convergence rate for the “bulged” 
decomposition is slightly higher than for the straight sided subdomains. This contrasts 
strongly with the observations of [36], which considered the approximation of second 
order problems. There, the presence of even slightly curved interfaces increased the error 
by orders of magnitude. 
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4 6 8 1012141618 
Polynomial Order 


Fig. 1 1. Convergence of the error for the three grids of Fig. 9. 

As our last example of the point source flow, we use the grid shown in Fig. 12. 
The geometry, a square with five circles cut out of its interior, was chosen to show that 
the method can be used to compute on a complex, multiply connected region. Twenty- 
tour subdomains were used to cover the computational domain. Up to seven subdomains 
share a common corner point without difficulty, because such points are not included in 
the discrete approximation. 
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Fig. 12. Grid for point source problem. 

The boundary conditions were chosen so that the exact steady solution was radial 
flow with the point source at the center of the middle circle of Fig. 12. The center cutout 
circle was specified as an inflow boundary, with the conditions chosen so that the Mach 
number of the incoming flow was M = 0.6. The boundary conditions along the remaining 
cutout circles were either inflow or outflow, depending on the direction of the normal 
velocity. The square outer boundary was an outflow boundary. For all inflow/outflow 
boundaries, the exact solution was used to provide the external flow values required by 
the Riemann solver. 

In Fig. 13, we plot the exact and computed Mach number contours for the solution 
of the point source flow. For the grid shown in Fig. 12, the contour lines of the exact 
solution, which are plotted with dashed lines, are coincident with the solid contour lines 
of the computed solution. 



Fig. 13. Solution of the point source flow for the goemetry shown in Fig. 12. The exact 
solution is plotted with dashed lines, the computed with solid lines. 

The approximation on the grid of Fig. 12 converges exponentially. Fig. 14 shows 
the maximum error in the density as a function of the polynomial order in each 
subdomain. We see that doubling the number of points per subdomain causes the error to 
decay by approximately two orders of magnitude. 

4.2 Subsonic Flow Over a Circular Bump in a Channel 

The second example is that of a Mach 0.3 subsonic flow over a circular bump in a 
channel. The geometry and grid with N = 9 is shown in Fig. 16. Wall boundaries were 
specified at the top and the bottom. At the left and right boundaries, the uniform flow 
free stream solution was specified as input to the Riemann solver. Initially, the free 
stream solution was specified everywhere, and then the boundary conditions were 
imposed. This problem does not have an exact analytic solution. However, since the 
incoming flow was chosen to be irrotational and isentropic, the entropy should be zero 
everywhere. The fact that this is not the case can be the result of the spatial 
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approximation and to the normal wave model used at the interfaces and boundaries for 
calculating the flux [32], 



468 10121416 

N 

Fig. 14. Convergence of the density error for the solution shown in Fig. 13. 


Solution contours of the Mach number for the grid shown in Fig. 15 are presented 

in Fig. 16. The wall pressure along the bottom, plotted as the pressure coefficient 
C p = (p- 1)/ 7 , is shown in Fig. 17. Finally, a convergence study of the entropy errors 

as a function of N is shown in Fig. 18. In that figure, we plot the maximum value of the 
quantity X = p/ p r - 1, which should be zero everywhere. We see that the error due to 
entropy generation converges exponentially fast. 
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Fig. 15. Geometry and grid for N - 9 for the flow over a circular bump in a channel. 





Fig. 18. Convergence of the entropy generated by the staggered grid approximation. 

4.3 Transonic Flow in a Converging-Diverging Nozzle 

As an example of a transonic problem, we compute the flow in an axisymmetric 
converging-diverging nozzle. We have chosen the nozzle used in the experimental 
investigation of Cuffel et al. [13], which was designed to show significant two 
dimensional effects. The nozzle consists of a converging section with half angle of 45° 
and a diverging section with half angle of 15°. The experimental tests were done in air 
with a stagnation temperature of 540 R and stagnation pressure of 70 psia. The nozzle 
geometry and the grid that were used in our computations are shown in Fig. 19. Note that 
we have varied the number of points per subdomain in this problem. 

To match the experimental conditions, we scaled the equations (31) and (32) by 
p = p* / p lot ,p = P* I P, ol * where the **’ represents the dimensional quantity. Under this 

scaling, the temperature and entropy are T = T‘ I T tot ,S lol = 0. The initial condition for 
the computation was the exact solution of the quasi-one-dimensional nozzle that has same 
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area as the two-dimensional nozzle. For the inflow condition at the left boundary, we 
specified that the tangential velocity be zero, the entropy be zero, and the temperature be 
unity. At the right boundary, the outflow is supersonic, so that no boundary condition is 
necessary. 



Figure 19. Grid for the 45° -15° converging-Diverging Nozzle 


Since not all of the external flow values are known at the left boundary, 
particularly the inflow velocity, it is not convenient to use (26) to impose the boundary 
condition. Instead, we use the following characteristic-like method that allows us to 
specify only the parameters that are known. The fact that the inflow condition sets v = 0 
means that the flow is essentially one dimensional. Then we can write a left-going 
Riemann invariant for the flow. In terms of the Mach number, M, and the sound speed, a , 
that invariant must satisfy 


f 2 ^ 


r 2 ^ 

l r-ij 

= R — n 

^ computed ^ comp 

M comp. j 

1 l J 


where the computed quantities represent values interpolated to the boundary by the 
reconstruction procedure. The boundary conditions fix the total temperature and the 
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entropy. With the scaling given above, this fixes the total sound speed at a tot = ^jy . Then 
the scaled sound speed and the Mach number are related by 


a = 


Vr 


1 + 



(52) 


Combining (51) and (52), an equation for the Mach number at the boundary can be 
written 






AT — 


7 - 1 , 


\ = R 


1 computed * 


(53) 


Equation (53) can be written as a quadratic equation in the Mach number and 
solved directly. Once the inflow Mach number is known, the sound speed can be 
computed using (52). From the Mach number, the sound speed, tangential velocity and 
the entropy, all remaining variables can be computed. From the full state on the 
boundary, the boundary flux can be evaluated. 

Results computed for the nozzle are shown in Figs. 20 - 23. Contours for the 
pressure are shown in Fig. 20. A comparison of the Mach contours and measured Mach 
number in the neighborhood of the nozzle throat is shown in Fig. 21. We see good 
agreement between the computed Mach contours and the measured values for Mach 
numbers up to about 1.6. We note that the discrepancies between the computed and 
measured Mach numbers are consistent with the discrepancies observed with the 
solutions of the inviscid flow solvers reported in [13]. Finally, in Figs. 22 and 23, we 
show a comparison between the computed and measured values of the pressure and Mach 
numbers along the upper wall of the nozzle. 
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Figure 21. Comparison of computed and measured Mach contours 
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Mach Number 



Fig. 22. Comparison of computed and measure wall pressure as a function 
of distance from the nozzle throat. 



Fig. 23. Comparison of computed and measured Mach numbers as a 
function of distance from the nozzle throat. 
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5 Concluding Remarks. 

We have presented a new, staggered-grid Chebyshev spectral multidomain 
method for the solution of inviscid compressible flow problems. The solutions are 
defined at the nodes of a Gauss quadrature rule, while the fluxes are evaluated at the 
nodes ol a Gauss-Lobatto rule. We have applied the method to one and two dimensional 
problems, but it should extend directly to three dimensions. 

The staggered grid multidomain method for compressible flow problems has 
many desirable features. These features include 

• Conservation. 

Mass, momentum and energy are conserved globally 

• Free-stream preservation. 

A uniform, steady flow stays uniform and steady even for complex 
subdomain shapes. 

• Temporal accuracy. 

• Geometric Flexibility. 

The domain need only be decomposable into quadrilaterals. 

• Programming simplicity. 

Corners of subdomains are not included as part of the approximation. 
Special cases do not need to be coded. Subdomains have at most four 
neighbors in two space dimensions and six in three space dimensions. 
Boundary conditions require no special comer treatments. 

While the new method is flexible, there remain a few limitations. We have 
assumed that the approximation is conforming. Thus, subdomains must meet at a point 
or along a full side. Along a side, the points at which the fluxes are computed must 
coincide with their neighbors. Also, the domain must be decomposable into 
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quadrilaterals. We expect that the first restriction can be eliminated by considering non- 
conforming approximations [31]. 

We have not considered the approximation of shocks in this paper. Methods for 
shock capturing with spectral methods that have been proposed to date include direct 
filtering [22], removal of the discontinuity plus filtering [6], Flux Corrected Transport 
methods [15], [16] and weak filtering with post-processing, [14], [17], Shock fitting 
should also be possible [23], [28]. A thorough study of the best options is beyond the 
scope of this paper. 

We have emphasized steady-state computations in this paper. However, the 
method is applicable to time dependent problems. Some one dimensional examples are 
included here and in [27]. 

The new method is a factor of two more expensive than methods that compute 
both the solution and the fluxes on the Lobatto grid. This is due to the extra interpolation 
from the Gauss points to the Gauss-Lobatto points. FFT techniques cannot be used to 
compute either the interpolation or differentiation operations with the new method. Thus, 
the approximation order within the subdomains must be kept low enough to where matrix 
multiplication is more efficient than FFTs, a value that varies from machine to machine. 
However, while the Lobatto grid methods require point operations at subdomain corners, 
the new method requires vector operations only. We believe that the flexibility and 
programming ease of the new method compensates for the additional work. 

The new method differs from the cell average method proposed in [35], [15], and 
[16]. The cell average method requires a pointwise procedure to be performed at 

subdomain corners that is not required by this method. The reconstruction and 
differentiation operations performed are different. In one space dimension the work 

required by the staggered grid method is twice that of the cell average method. However, 
the reconstruction procedure for the cell average method in two space dimensions 
requires two matrix multiplication operations per line in each direction, as does the 
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staggered grid method, so the work is equivalent in two space dimensions. In three space 
dimensions, the staggered grid method will require only two thirds the work of the cell 
average method, making it more efficient in that case. 
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